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Based on a theoretical argument and Monte Carlo simulations of a Ginzburg-Landau 
model derived microscopically, it is argued that, in type-II superconductors where both 
the paramagnetic and orbital depairings are important, a strong first order transition 
(FOT) at H C 2 expected in the mean field (MF) approximation never occurs in real 
systems and changes due to the fluctuation into a crossover. The present result explains 
why a nearly discontinuous crossover at H C 2 with no intrinsic hysteresis is observed 
only in a clean superconducting material with a singlet pairing and a high condensation 
energy such as CeCoIn 5 . 
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The phase diagram of a type-II superconductor in nonzero magnetic fields {H ^ 0) has been 
studied in relation to high T c cuprates in fields of the tesla range and beyond the mean field (MF) 
approximation completely neglecting the superconducting (SC) fluctuation effect. 1, 2 ) Even in a 
material with a spin-singlet electron-pairing, the spin (i.e., Pauli paramagnetic) depairing effect 
can be neglected in relatively low fields, and the transition at H C 2(T)-line in the MF approximation 
for this case is of second order. On the other hand, the phase diagram in cases where both the 
Pauli paramagnetic and orbital depairings are important were examined so far only within the MF 
approximation. 3 ' 4 ) Since the MF transition in this case has been expected to become first order at 
low enough temperatures, 3, 4 ) it is valuable to clarify whether this MF scenario is also, as well as 
in the low field case, 1,2 ) drastically changed by including the SC fluctuation. 

In this paper, we first examine the MF approximation of a quasi two-dimensional (2D) su- 
perconductor under a field perpendicular to the SC layers {H || c) by treating the orbital- and 
spin-depairing effects on an equal footing and obtain a first order transition in the MF approxi- 
mation (MF-FOT) at H C 2(T) and in high fields. Then, based on a physical picture on the high H 
fluctuation and our Monte Carlo simulation results on the derived Ginzburg-Landau (GL) model, 
we argue that any nonvanishing SC fluctuation, as well as in the low H case, 1,2 ) changes this MF- 
FOT into a crossover and point out a consistency between the present theory and experimental 
results, such as the apparently discontinuous behavior 5-8 ) in high fields with no intrinsic hysteresis 
seen in CeColns with unusually large condensation energy and the absence of such a sharp behavior 

* ada@scphys.kyoto-u.ac.jp 



1 



2 Hiroto ADACHI, Shigeru KoiKEGAMl,and Ryusuke Ikeda 

in other materials 9-11 ) showing large paramagnetic effects and having lower condensation energies. 

Let us first briefly explain our microscopic derivation of a Ginzburg-Landau (GL) functional 
in the case with spin depairing effect by starting with a quasi 2D BCS hamiltonian. An attractive 
interaction leading to the d x i _ y i -pairing state will be assumed. For H \\ c case, the following GL 
model will be derived microscopically for a layered system: 



T = N(0) J d 2 r ± 



q z N 



+ £(^|a?W + ^|a<°W) 



(i) 



where Aj N \r±) = N^ 2 J2q z Aq^\r±)e iqzS: > denotes the pair-field (SC order parameter) on the 
j-th SC layer projected onto the Nth Landau level (LL), N(0) is a 2D density of state at the Fermi 
surface, -/V"i ayer is the number of SC layers, and s is the layer spacing. Below, the higher LL modes 
with N > 2 will be neglected. In eq.(l), the 4th and 6th order terms were represented within the 
lowest LL approximation and in a spatially local form. The neglect of this non-locality is justified 
because the energy scale associated with the MF transition (i.e., the condensation energy) is much 
larger than that on a phenomenon arising from this non-locality such as a structural transition 
between the rhombic and square vortex lattices. Then, assuming, for brevity, a cylindrical Fermi 
surface and treating effects of the magnetic field quasi-classically, the coefficient a^iq 2 ) of eq.(l) 
in pure case (r _1 = 0) becomes 

a N {q 2 z ) = HT/T cQ )-27rTj2(D N (2e)-l/e), (2) 

£>0 

D N (2e) = 2 jH dpe- 2ep -(^ C N (2{^-) 2 ) 

xcos(2n Hp)J (2Jsm{^)py (3) 

where e = 2ir(n + 1/2)T is a Matsubara frequency, T c q the transition temperature in H = 0, th = 
th/vf, vf the Fermi velocity, rn = y / (f)o/2irH the magnetic length, J the interlayer hopping, ^qH 
the Zeeman energy, and Cn and Jn are the Nth order Laguerre and Bessel functions, respectively. 
In the impure case with a finite elastic scattering rate r _1 , the impurity-ladder vertex correction 
can be neglected in non s-wave pairing cases if th *C t. 3 ' 12 ) Then, the corresponding expression 
to eq.(2) is given by replacing D^{2e) by Dn(2e + r -1 ). 

The corresponding calculation of the coefficients V4 and Vq is more involved and requires mul- 
tiple numerical integrals. The orbital depairing always introduces the field-induced vortices, which 
may change drastically the phase diagram in H / 0, and hence, will be included nonperturbatively 
in contrast to ref. 4 ) Further, the paramagnetic effect becoming important upon cooling is sensi- 
tive to the impurity scattering. Hence, both the two depairing effects and the impurity scattering 
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Fig. 1. (a) Mean field H-T phase diagram obtained from computations in high H and low T region. The magnetic 
field is measured in the reduced variable h = ff/ff°| blt . (b) T-dependences of V4 = Va/t^ and V& = Vq/t% 
computed on H C 2{T). In both (a) and (b), (2-ktTco)^ 1 = 0.025 was commonly used. See the text for other 
details. 

should be treated on an equal footing. For this reason, we have examined the coefficients V4 and Ve 
numerically by, as well as in eq.(2), applying the parameter (p in eq.(3)) integral representation and 
neglecting the vertex corrections. Details of their derivation will be described elsewhere. 12 ) Conse- 
quently, as shown in Fig.l, the coefficient V4 (Vq) defined in N = becomes negative (positive) in 
H > H*(T). 

In Fig.l (a), we give a MF phase diagram in high H following from calculations of the GL 
coefficients mentioned above in a slightly impure {(2ittT c q)~ 1 = 0.025) case. There, the param- 
agnetic parameter fioH°2 blt /2ttT c q corresponding to the Maki parameter is chosen as 0.8, where 
H™ hit = O.560o/(27r£g) is the 2D orbital limiting field in pure case, and £0 = vf/2ttT c q. The 
dashed curves express the lines defined by ao(0) = and a±(0) = 0, and, in a range where the 
MF transition is of second order, the higher field portion of these two curves plays the role of the 
H C 2(T)-\me. Above the curve denoted by the plus symbols (i.e., H > H*(T)), however, the coeffi- 
cient V4 in N = is negative, making the MF transition in JV = a discontinuous one (i.e., FOT) 
which occurs on the portion of H&iT) denoted by the open circles in T ncxt < T < T*. On the other 
hand, the MF transition in N = 1 is found to remain continuous, 12 ) although the N = 1 region 13 ) 
with ao(0) > ai(0) is quite narrow. We have also estimated a structural transition curve //fflo(T') 
between the ordinary straight vortex solid and a FFLO-like helical vortex solid with a modulation 
A^- ~ e iqzS ^ along H as the line where da^ql) / dq^\ qz= o = 0. Since 9 2 ao((/ 2 )/9((? 2 ) 2 |g z =o > 0, 
and the spin depairing inducing this transition is enhanced with increasing H and decreasing T, a 
second order structural transition at H-fflo{T) decreasing upon cooling is expected at least near 
H C 2 to occur, although the range in which this N = state is realized seems to be quite narrow in 
this H || c case. Actually, in contrast to T*, the i^FFLo(^) curve and the N = 1 vortex state region 
are very sensitive to the impurity strength, and we find that, for {2ttT c qt)^ 1 = 0.05, T* ~ 0.23T c o, 
while the H-pYhoiT) and the N = 1 region are not present any longer. The main result T* > Tfflo 
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in our MF calculation is satisfied, at least in H \\ c, irrespective of the impurity strength 12 ) and 
should be compared with results in ref. 4 ) 

Now, let us turn to discussing the genuine high H phase diagram of superconductors with 
paramagnetic depairing. Hereafter, we will focus on the region H*(T) < H < Hp-phoiT) in 
Fig. 1(a) and use the following GL functional on the basis of the MF results shown above: 

f = £|dV ± {a|tf/ + 7 |tf.;-.p J+1 | 2 

+ (4) 

Here, 7 > 0, a ~ ao(H — Hq)/Hq with a T-dependent constant ao, and Ho is the field value 
satisfying ao(0) = 0. Further, the length and the pair-field were rescaled in the manner r±/rjj — ► r± 
and (T/ Vq r 2 ^iV(0))~ 1 / 6 A^ ' ) — ► iffj. Consequently, the 2D fluctuation strength is measured by 
\f3(T)\~ l = \V4\~ 1 (VqT / N (Oyjj) 1 / 3 . Any term expressing a vortex pinning effect induced by the 
electronic impurity scattering was neglected. 

To explain why the MF-FOT should not occur in real systems with any nonvanishing fluctua- 
tion, it is convenient to start from a description deep in the ordered phase. First, since an inclusion 
of the orbital depairing requires the presence of field-induced vortices, a low energy excitation in 
the ordered phase is inevitably an elastic mode of a vortex solid. It is quite clear that the form of 
the elastic energy and the relation between the phase fluctuation and the vortex displacement 14 ) 
are unaffected by a difference in the form of nonlinear terms of the GL model. As in the case with 
no paramagnetic depairing, 14 ) the phase fluctuation is marginally relevant even in 3D case, and the 
main rigidity of the phase fluctuation is the shear modulus of the underlying vortex solid. Hence, 
if the melting of the vortex solid occurs at H m below H C 2, the vortex liquid just above H m with 
vanishing shear modulus has no long range orders and hence, should be continuously connected 
with the normal phase above H&P^ In this case, the MF-FOT at H C 2 should not occur. Further, 
to understand this from another point of view, let us note that the quasi 2D SC order parameter 
in the lowest LL has the form 16 ) 

N 3 -l 

z) = A{z)e- y2 l {2r ^ H (u - Ui (z)), (5) 
i=0 

where u = x + iy, z = js, Ui(z) is the complex coordinate perpendicular to H of the i-ih vortex, 
and a Landau gauge was assumed for the external gauge field. Since the vortex positions are highly 
disordered above H m , the fluctuation effect above H m is essentially described only by the amplitude 
A(z). However, A(z) depends only on z irrespective of 'the form of nonlinear terms in the GL model. 
That is, since the amplitude fluctuation itself has a reduced dimensionality and is ID-like in 3D 
systems 2 ' 15 ) even in the present case, the amplitude fluctuation plays a role of pushing the vortex 
lattice freezing field down to a lower field, and the MF-FOT should not be realized in real 3D 
systems. 
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Fig. 2. Numerical data (symbols) of (l^t 2 ) in the case composed of four layers. The corresponding results (curves) in 
the Hartree approximation 15 - 1 are also given. The cross symbols and solid curve (plus symbols and dotted curve) 
are results in |/3| = 2.0 (3.0). In \f3\ = 3.0, the H C 2(T) corresponds to a = 1.7 and lies above H m - The used 
parameter values are given in the text on the final page. 

Bearing such a theoretical prospect in mind, we have performed Monte Carlo simulations on 
the model of eq. (4) , and its 2D version using material parameter values useful for comparisons with 
observations in CeCoIns. Our simulation method follows that of ref. 17 ) under the lateral boundary 
condition L x /L y = 2N x /y/3N y with (N x ,N y ) = (6,6) and a periodic boundary condition across 
the layers. Hereafter, the key quantities in the present case within N = are — (|^| 2 ) and the 
structure factor defined 15 ) as the correlation function of |^(^j_)| 2 - The former is equivalent to the 
magnetization or the fluctuation entropy, while the latter is a measure of the vortex positional 
order. 




Fig. 3. Structure factor defined from the correlation function of |\&(r)| 2 for |/3| = 3.0 at (a)a = 1.65,(b) a — 1.25. 
Computed results of a dependences of (\^\ 2 ) in a layered system consisting of four layers are 
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shown in Fig.2 for \f3\ = 2.0 and 3.0. Clearly, the nearly discontinuous behavior at H C 2 is smeared 
out as the fluctuation is enhanced (i.e., \(3\ is reduced). For comparison, the corresponding results in 
the Hartree approximation 15 ) are also drawn as the solid and dashed curves. Note that a change of 
(l^l 2 ) at H C 2 is extremely steep even in this approximation where no transition occurs. 15 ) Clearly, 
the presence of a genuine FOT cannot be argued merely through a steep growth of (|^| 2 ) in a 
real system with weak SC fluctuation. Next, let us examine whether the melting transition point 
coincides or not with the MF transition field H C 2(T). The melting transition is widely believed to 
be a weak FOT, and this should be found 17 ) in Monte Carlo simulations as a tiny discontinuity in 
thermodynamic quantities. Unfortunately, due primarily to numerical difficulties, our simulation 
is restricted to too small systems to observe such a discontinuity. For our purpose, however, it 
is sufficient to estimate the melting temperature by finding where the Bragg peaks of the vortex 
lattice disappear. Fig.3 shows snapshots of the structure factor on |^| 2 when \(3\ = 3.0. No vortex 
positional ordering is seen at H C 2. As indicated in Fig.2, nearly sharp Bragg peaks appear at H m 
below H C 2, while most of the entropy have been lost rather near H C 2 above it. The two field (or 
temperature) scales, one characterizing the steep growth of (l^l 2 ) and another corresponding to 
the sudden growth of vortex positional ordering, are clearly distinguished. 

Quite similar results to those in Fig.2 and 3 were also obtained in 2D case where the system 
size dependence of H m is more easily examined. As expected from a commensurability between a 
MF solution and the lateral periodic boundary condition, it was verified through the data in system 
sizes (6,4) and (8,6) with L x /L y = \/3N x /2N y that H m is lowered with increasing the system size. 
This means that the relation H m < H C 2 found above holds in simulations for larger sizes. 

Recently, a large apparent discontinuity was observed in transport and thermodynamic data 6-8 ) 
in a heavy fermion material CeCoIns in high fields, where a strong paramagnetic depairing is ex- 
pected even in H || c, 6 ' 18 ) and was identified with the MF-FOT. A hysteresis seen in the magneti- 
zation data was interpreted as being due to a genuine FOT at H C 2. However, the specific heat in 
H || c 8 ) and resistivity data 7 ) have shown no measurable hysteresis in contrast to observations 19 ) 
in systems with no orbital depairing, although a stronger genuine FOT should be accompanied 
by a larger hysteresis reflecting a larger condensation energy. Further, a magnetic hysteresis in 
type-II superconductors is usually a consequence of vortex pinning. 20 ) Actually, a magnetization 
measurement in CeCoi-zRhjJns has shown 9 ) that the discontinuity of magnetization disappears 
with increasing the Rh-substitution, while the hysteresis remains unaffected, implying that the 
observed 6 ' 7 ) magnetic hysteresis is not a direct consequence of an intrinsic SC transition. Further, 
the Rh-substitution has resulted in a decrease of condensation energy 9 ) which is usually equivalent 
to an enhancement of SC fluctuation. Together with the absence of a nearly discontinuous behav- 
ior in other systems 10 ' 11 ) suggestive of a large spin depairing and with strong fluctuation, these 
observations in CeCoIn 5 are consistent with the absence of a genuine FOT stressed in this paper. 
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Bearing these observations in Ce-compounds in mind, we have chosen parameter values in Fig.2 
and 3 as follows. The iV(0)-value was taken from the heat capacity data in H = 0. 18 ) The Maki 
parameter value used in Fig.l leading to T*/T c q ~ 0.37 was chosen again because the apparent 
FOT behavior is observed in T < 0.3T c q. 8 ^ Further, we have assumed, for brevity, T/T c q = 0.1 and 
7 = 0.25 (see eq.(4)). Then, the -value is located in the range between 2.0 and 3.0 used in our 
computation, and we have a (oc (N(0)T c0 /T) 2 / 3 ) ~ 50. This a value suggests that, on a linear 
scale of the magnetic field of the order of H^iT = 0), even the broader simulation curve in \(3\ = 2 
in Fig.2 cannot be distinguished from a strictly discontinuous jump. Further, according to our MF 
calculation the (apparent) magnetization jump value at T = 50(mK) is roughly estimated as 20(G) 
which is comparable with the measured one (~ 30(G)). 6 ) 

In conclusion, we have given convincing results implying that a MF-FOT at H c2 in supercon- 
ductors with field- induced vortices (i.e., orbital depairing) should not occur as a genuine FOT in 
real systems with any nonvanishing fluctuation effect. Typically, this situation is realized in mate- 
rials with a singlet pairing and with an enhanced paramagnetic depairing. Assuming a quasi 2D 
material with the d x i_ y i -pairing state in H || c, we have performed numerical simulations under a 
condition presumably compatible with high field CeCoIns data suggesting a strong discontinuous 
transition at H C 2. The absence of a genuine FOT at H C 2 is consistent with the absence of hysteresis 
accompanying this discontinuous behavior in CeCoIn 5 and unrelated to the vortex pinning. The 
obtained MF phase diagrams seem to be roughly consistent with the observations in CeCoIn 5 , 
including recent data 21 ) suggestive of a structural transition possibly to an FFLO-like vortex solid, 
if taking account of a nonvanishing impurity concentration. The issues in H _L c will be considered 
elsewhere. 

After writing an earlier version of the manuscript, additional simulations at lower T/T c q (< 0.1) 
have been performed, and a small hysteresis in the almost discontinuous behaviors was found to 
occur in this case. However, since the GL coefficients (i.e., V4 and Vq in eq.(l)) remain unchanged, 
this hysteresis is not due to a genuine FOT but a phenomenon appearing on measurable Monte 
Carlo steps (i.e. time scales) as a reflection of incomplete relaxation in simulations 12 ) which are 
inevitably in non-equilibrium. The tiny hysteresis appeared at lower T but above Tfflo in CeCoIns 
only in H _L c 21 ) is likely to have a similar origin to this numerical one practically occurring in 
cases with sufficiently weak fluctuations. 

We thank T. Sakakibara, Y. Matsuda, K.Izawa, T. Tayama, R. Movshovich, and K. Machida 
for stimulative discussions. Part of the numerical computation in this work has been carried out 
at the Supercomputer Center, Institute for Solid State Physics, University of Tokyo, and at the 
Yukawa Institute Computer Facility in Kyoto University. 



1) E. H. Brandt: Rep. Prog. Phys. 58 (1995) 1465. 



8 Hiroto ADACHI, Shigeru KoiKEGAMl,and Ryusuke Ikeda 

2) R. Ikeda: J. Phys. Soc. Jpn. 70 (2001) 219. 

3) V. R Mineev: Phil. Mag. B 80 (2000) 307. 

4) M. Houzet and A. Buzdin: Phys. Rev. B 63 (2001) 184521. 

5) K. Izawa et al: Phys. Rev. Lett. 87 (2001) 057002. 

6) T. Tayama et al: Phys. Rev. B 65 (2002) 180504. 

7) T. P. Murphy et al: Phys. Rev. B 65 (2002) 100514(R). 

8) A. Bianchi et al: Phys. Rev. Lett. 89 (2002) 137002. 

9) A. Harita et al: in Autumn Meeting of the Physical Society of Japan (2002). 

10) T. Shibauchi et al: Phys. Rev. Lett. 86 (2001) 5763. 

11) T. J. Lee et al: Phys. Rev. Lett. 78 (1997) 3555. 

12) H. Adachi and R. Ikeda: in preparation. 

13) H. Shimahara and D. Rainer: J. Phys. Soc. Jpn. 66 (1997) 3591. 

14) R. Ikeda et al: J. Phys. Soc. Jpn. 61 (1992) 254. 

15) R. Ikeda et al: J. Phys. Soc. Jpn. 59 (1990) 1397. 

16) Z. Tesanovic and L. Xing: Phys. Rev. Lett. 67 (1992) 2729. 

17) Y. Kato and N. Nagaosa: Phys. Rev. B 48 (1993) 7383. 

18) S. Ikeda et al: J. Phys. Soc. Jpn. 70 (2001) 2248. 

19) T. Suzuki et al: J. Phys. Soc. Jpn. 69 (2000) 1462. 

20) A. A. Gapud et al: Phys. Rev. B 67 (2003) 104516. 

21) A. Bianchi et al: cond- mat/0304420. 



